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Abstract 

A particular, two-dimensional, tiling model, composed by the so-called 
Wang tiles has been studied at finite temperature by Monte Carlo nu- 
merical simulations. In absence of any thermal bath the Wang tiles give 
the opportunity of building a very large number of non-periodic tilings. 
We can construct a local Hamiltonian such that only perfectly matched 
tilings are ground states with zero energy. This Hamiltonian has a very 
large degeneracy. The thermodynamic behaviour of such a system seems 
to show a continuous phase transition at non zero temperature. An order 
parameter with non-trivial features is proposed. Under the critical tem- 
perature the model exhibits aging properties. The fluctuation-dissipation 
theorem is violated. 

1 Introduction to the model 

In this work we propose a model that shows a non-trivial thermodynamic be- 
haviour due to its particular geometric structure. 

This two-dimensional model is built using square tiles, called Wang's tiles, 
on a square lattice. The edges of these tiles can be of six different 'colours'. 
But of all the possible types (6^) that can be created changing the colours just 
a particular group of sixteen, found by Ammann |^] , is taken into account (see 
figure p. 

This is one of the minimal set of Wang tiles such that the corresponding 
tilings exist and are non-periodic tilings. A tiling is a configuration of tiles 
placed edge-to-edge on the plane, where all contiguous edges have the same 
colour. If at least one tiling is allowed and none of the possible tilings shows a 
periodic pattern, then the set of tiles composing them is called aperiodic. Other 
aperiodic sets of tiles are often used as models for quasi-crystal materials 
but this is not the present case. 

The Wang tiles were the first aperiodic tiles to be discovered, in 1966 by 
Berger p|. They were initially important, and they still are, because of their 
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Figure 1: The sixteen Wang tiles. The six different types of edges, or 'colours' are 
indicated here by numbers (0, 1, 2, 3, 4, 5). The tiles shown are the minimum set 
allowing aperiodic tiling 



use in problems of mathematical logic [yQ. The use that we make of them 
in this work is, anyway, far away from that point of view. Indeed we look 
at the behaviour of a system built by Wang tiles in a thermal bath, defining 
thermodynamic observables on these tilings. The main stimulus to study this 
model has been the high degeneracy of the perfectly matched configurations and 
their non-trivial aperiodic structure. 

Putting the system in a thermal bath allows for translations of the tiles 
also in positions where there is no matching between edges of neighbours tiles, 
forming in this way a unmatched-tiling. We can assign an energy to each con- 
figuration which is equal to the number of links such that the colours of the 
facing edges are different. The energy of the exactly matched configurations, 
that from now on we will call ground states, is then zero. 

If we label the type of tile in the position given by the coordinates {x, y) in 
the plane by T^^y and the type of its four edges respectively towards South, East, 
North and West by T^^y , T^^y , Tx^y^ and li^^ we can write an Hamiltonian for 
this system in the following way: 

{x,y) {x,y) 

L is the number of tiles in the lattice in one direction. The 5{z) is the delta of 
Kronecker. 



2 Equilibrium analysis and critical behaviour 

The numerical simulations done in order to study the equilibrium characteris- 
tics have been performed using the parallel tempering algorithm p]. We have 
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Figure 2: Left: specific iieats in the wiiole probed region. Rigiit: specific iieats 
around tlie critical region. Linear sizes from 8 to 32 are plotted. The crossing point 
shift with the size of the system following a FSS. We can observe both the raising of 
the slope of the curves around the crossing point and a growing peak at T = 0.42 



chosen open boundary conditions on the two dimensional lattice. Indeed using 
periodic boundary conditions the fact that periodic tilings do not exist would 
have implied that the energy of the ground state would have been different from 
zero. With open boundary conditions the energy of the ground state is, by 
definition of (|l|), equal to zero. It is equivalent to say that all the tiles are edge 
by edge perfectly matched, forming some non-periodic structure. 



2.1 Phase transition 

At different temperatures, for every system we have computed the energy and 
the specific heath. Numerical simulations have been brought about on square 
lattices of different size: from a linear size of 8 to one of 32. Every equilibrium 
simulation has been carried out for a number of MC steps going from 10 millions 
to 100 millions, relatively to the size of the lattice. A range of temperature 
between and 5 has been observed, then concentrating in small intervals in 
what comes out to be the critical region, around T = 0.4. We controlled that at 
zero temperature the energy goes to zero and that the specific heat computed 
from the derivative of the energy coincides with the one computed from the 
energy fluctuations. These are checks of a correct thermalization. 

The specific heat presents a smeared but clear change around the tempera- 



ture 0.4 (see figure 2.1), from a lower value for lesser temperatures to a higher 
value for greater. Increasing the size of the system we can observe that the 
crossing points between two curves of different size moves towards right and 
that the slope of the Cl (T) line increases in this interval of temperatures as L 
increases, but we also see that a peak grows around T = 0.42. 

Performing Finite Size Scaling (FSS) analysis of the crossing points of the 
specific heat curves we find that their abscissa tends to a Tcrossioo) = 0.398 ± 



3 



0.007 as L — + oo, with an exponent v = 1.6 ± 0.5 and that the behaviour of the 
slope with increasing size is compatible with a diverging fit at that temperature. 
In this case we would have a jump at Tc = Tcross(cx)), which corresponds to a 
critical exponent a equal to zero. 

But at the same time the fit of the peak height is consistent with a divergence 
at Tc = 0.42 ± 0.01. Knowing that at criticality C - |T - Tc|~" and ^ ~ 
|T — Tc\~^ ^ L we get C ^ L"/". Our data are both consistent with a power- 
law divergence {a/iy — 0.35 ± 0.01), and with a logarithmic one (a = 0). A 
similar value of the critical temperature was found, in this last hypothesis, by 
Janowsky and Koch ^ . 

In any case there is evidence for a second order phase transition. 

2.2 Order parameter 

Wang's tiles satisfy a particular property: if all the tiles of the two contiguous 
west and south sides (or equally of the north and east sides) of the square lattice 
are fixed, then at most one aperiodic tiling can be formed. This is due to the 
fact that each tile of the Wang set has different combinations of west-south 
(or north-east) colours. Then at most only one of them can be put, in the 
bottom-left (either top-right) corner. The same is true for the tiles to be put 
in the two corners created by putting the first tile. If any, there will be only 
one combination available also for them. Of course this is valid at T = in our 
system. 

This condition implies that ground state degeneracy can not increase faster 
than 16^^^^ and therefore the entropy density is zero at zero temperature: 
sl(0) a/L, where a < 81og2. Since just a few north-south or west-east com- 
bination are allowed the actual constant in the entropy is smaller than 8 log 2 
(a stricter upper bound of a = log 12 = 3.585 log 2 can be easily obtained). 

In order to identify an order parameter we have to break this degeneracy. 

We have then simulated the parallel evolution of two copies of the system, 
with the following procedure: one system reaches the equilibrium, then a copy 
of the equilibrium configuration is done and the evolution of this second copy is 
performed. One foundamental constraint is set: the boundary tiles on the south 
and east sides of the lattice stay unchanged in the dynamics towards equilibrium 
of the second copy. 

For r > the equilibrium states are no more the exactly matched tilings. 
Because of the thermal noise some couple of neighbours sides of the Wang's tiles 
can be unmatched. Furthermore there won't be a unique tiling minimizing the 
energy. The second copy can then evolve to a different configuration. 

We are interested in looking what happens when the temperature is increased 
from below to above the critical temperature, and how the behaviour is sensitive 
to the size of the system. With this aim we introduce an overlap that depends 
on the distance I of the tiles from the two contiguous boundary sides that are 
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fixed once that the first copy has reached equilibrium: 



9(0 



2{L-l) + l 



(1) 
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y>l^x=l x>lyy=l 



(2) 



Where the average (...) has been performed over different realizations of the 
configurations, ((...)) is the time average at equilibrium and I = 1, ...,L. 

The overlap q{l) is computed along the diagonal: starting from the vertex 
shared by the two fixed contiguous sides, q{l), and ending at the opposite vertex 
of the lattice, q{L). In order to gain more statistics the overlap is built averaging 
over the 2{L — l) + l elements in the row of ordinate I with the coordinate x > I 
and in the column of abscissa I with y > I and assigning this average value to 
the 'diagonal' function q{l). 
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Figure 3: q{l) profiles at different temperatures and difTeront sizes. T goes from 0.2 
above to 5.0 at the bottom. The behaviours for L = 16, 20, 24, 32 are shown. For 
L = 16 also the behaviour at T = 10~® is plotted. 



5 



The statistic sample is formed repeating the simulation several times with 
different initial configurations, in order to obtain different equilibrium config- 
urations (we have always checked if the same equilibrium configuration could 
possibly appear more than once starting from different initial configurations to 
avoid annoying biases of the statistic sample, but this never happened). Every 
size has been simulated for at least 100 different initial configurations. 

At zero temperature the overlap is always one. 

At higher temperature but below the phase transition it goes to a constant 
lesser than one (far enough from the boundary L, where finite size effects are 
overwhelming) . 

Above the transition point q{l) decays very rapidly, compatible with an expo- 
nential, to the lowest possible value, corresponding to completely uncorrelated 
copies {hot or disordered phase). Since the tiles are of sixteen different types 
and the probability distribution of the tile-types is uniform, this lowest value is 
equal to 1/16. 

In figure ^ we show the behaviour of q{l) at different temperatures, both 
above and below the critical one, for different sizes. From this figures we can 
see that the q{l) seems to approach some plateau for / L/2 at temperatures 
below T ~ 0.42. This is more evident as we go to bigger sizes. 

We can compute a q{T) for every size, averaging over the plateau values of 
g(/). The plateau is every time chosen taking a small window around L/2 and 
then enlarging it as long as the plateau value stays constant. As I grows further 
finite size effects destroy this plateau. We observe that this parameter shows a 
change approaching Tc. when T < Tc, it increases from the value of 1/16 that 
it has at high temperature to the value of 1 that it reaches at zero temperature. 
The profiles of qL{T) are plotted in figures ^. 




Figure 4: qL{T) profiles for L = 16, 20, 24 and 32. qL{T) is tlie plateau value computed 
respectively on intervals of 3, 4, 6 and 10 tiles 
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Figure 5: Probability distribution of the overlap's values for three different temper- 
ature: T = 0.45 is already in the disordered phase, while T — 0.35 and T — 0.2 are 
under the critical temperature. The behaviour is slightly dependent on temperature. 
For T — > the distribution tends to a (5 function on the value 1. 



Figure 6: Different profiles for different realizations of the initial conditions at T = 
0.35 for a system of linear size 32. 



The probability distribution of the overlap's values at the plateau has a 
non-trivial shape as soon as T < Tc (figure ||). This is not so astonishing since 
fixing the boundary conditions we have broken the degeneracy of the equilibrium 
states. (We stress that due do the overlap's definition the distribution doesn't 
show any symmetry q —>■ —q). The form of the P{q) averaged over different 
realizations of the tile's configurations at two continuous edges of the lattice is 
not strongly dependent from the size of the samples, at least for the simulated 
cases. 

The probability distributions of the overlap for different realizations of the 
equilibrium boundary conditions in the bottom and left edges of the tiling are 
represented in figure ^. 
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Figure 7: 'directions' of tlie random field fe'^*"^ in tlie space of tile types, represented 
in the complex plane. 



3 OfF-equilibrium analysis 

In order to study systems with very slow relaxation times, much greater than the 
experimental times, it is very important to look at the ofF-equilibrium dynam- 
ics. This approach describes glasses, spin-glasses and, in general, any system 
with many equilibrium states divided by high free energy barriers, in a more 
realistic way in comparison with an equilibrium approach, since equilibrium is, 
in practice, never reached during experiments. 

To keep the sample out of equilibrium during the run of a numerical sim- 
ulation we need the typical distance ^(i) over which the system has reached 
equilibrium at a certain time to be always smaller than the linear size L of the 
lattice. This distance is also called dynamical correlation distance. Practically, 
in order to be sure of avoiding thermalization, we choose sizes bigger than those 
used in the static analysis. The greater is the size the longer is the time needed 
to the correlation distance to reach the size of the system. More over we used in 
this case a standard Monte Carlo algorithm instead of the parallel tempering, so 
that the thermalization times increased sensitively under the transition point. 

In order to study the response of the system to an external field we can 
embed the system into a perturbative field. Namely a field directed in one of 
the sixteen possible directions in the 'tile-types space', where a particular tile of 
the set of sixteen can be favoured (positive field) or disfavoured (negative field). 
This can be a uniform or non-uniform field. To avoid any preference towards 
one particular type of tile we have chosen a non-uniform random field, whose 
value at each site is independent from the other sites. 

Thus we add to the Hamiltonian (1^) the perturbative term: 



where h^.y = h^xy'^^xyi'' ^he random external field pointing along one of the 
tile- types or opposite to it. 

h*^'^'^ has a uniform distribution of the sixteen possible choices of tile-type. 




(3) 
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and h%^^" gives the magnitude of the field and its sign (it can be randomly 
positive or negative, always according to a uniform distribution). 
The probability distribution can then be written as: 

1 1 

nh.,y) = Y^T.^ (^^T - e'"^) X 2 ^^^^^ - + '^('^^'r + ^°)) (4) 

n 

with hx.y = and h'^ y = h^- n = 0, 15 gives the 'direction' of the field in a 
representation in the complex plane (see figure ^ 

Once the system is cooled from high temperature, it is left evolving until a 
certain time f^, usually called waiting time. Kt t = tw the field is switched on 
and we begin recording the values of the temporal correlation function C(t, tw) 
and of the integrated response function m{t,tw)- 

For our model they are defined as follows: 

C{t.tu,) ^^Y.^{T.'y(t)^T,,y{t^)), (5) 



1 ^(/4T(^..) s{T,jt)-h'yr{t^))) 

m{t,t^;h)^-^J2- J, 



L2 



and the susceptibility is 



, ^ m[t,t^;h) 
X(t, t^) = hm . (7) 

ho— IIq 

For numerical computation it becomes: 

X(M.)-^^^^^. (8) 
ho 

Here ((...)) is the average over different dynamical processes and (...) the average 
over the random realizations of the perturbative external field. In a system at 
equilibrium the Fluctuation-Dissipation Theorem (FDT) holds: 

X{t-t^)=^—^^. (9) 

where we have made explicit use of the fact that the correlation function is 
defined in such a way that C{tw,tu,) = C(0) = 1. 

Out of equilibrium this theorem is no more valid. Anyway a generalization 
is possible 0, at least in the early times of the dynamics. This generalization is 
made introducing a multiplicative factor X(t,t„) depending on two times such 
that : 

X(t,i^) = ^^|^(l-C(t,i^)) (10) 

In a certain regime, called aging regime, X{t, t^) < 1 and the FDT is violated. 
An important assumption is that this modified coefficient X(t,tw)/T depends 
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on t and only through C{t, t^j). Its inverse is also called effective temperature 
Tf. = T/X [C{t,tuj)] since the system in this regime, for a given time-scale, 
seems to behave like a system in equilibrium at a temperature different from 
the heat-bath temperature. 

In terms of susceptibility and correlation functions we can write a general 
functional dependence: 

x{t,t^)^^S[Cit,t^)] (11) 

The <S'[C] here defined would be 1 — C if we were at equilibrium. 

From our probe we find that we first have a regime where the relation is 
linear with a coefficient equal to the heat-bath temperature. This early regime 
is sometimes called stationary, since the observables computed at very short 
times (compared with t^) do not depend on the age of the system. During this 
regime the system goes fast towards a local minimum. 
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Figure 8: L=128, h=0.2, S[C] at different temperatures. T — 0.45 is above the phase 
transition: in this case there is no regime in which the fluctuation-dissipation theorem 
does not hold. The time evolution of the system for t > has to be read going from 
right (C(i„,t„) = 1) to left. 



After this first time the S[C] bends and the coefficient of X{t,tw) /T is no 
more the inverse heat-bath temperature. X{t,tio) is now less than one, as if 
the system would be at an effective temperature bigger than the one of the 
heat-bath. It also seems to change continuously as the system evolves, until it 
reaches zero (figures ^ - p^ . 

The value of the autocorrelation function at which S[C] leaves the 1 — C 
line, equal to the average overlap value q in the static, increases continuously 
with temperature, as we already observed in the static analysis. 

The S[C] shows no strong dependence from the magnitude of the field, at 
least for the values that we have used to perturbate the system. Furthermore 
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Figure 9: Left: T=0.2, h=0.2, S[C] for different sizes. The dependence from the size 
is little. Right: L=256, T=0.2, S[C] for different values of h. 



the system shows only a slight dependence from the size (measures on L = 
64,128,256) (fig. ||). 

From the figures |l^ we can see that S[C] moves towards some asymptotic 
line increasing tw 

The initial difference between 1 — C and S[C] is due to lack of statistics: 
if Nfr is the number of field realization performed in the simulation there is a 
difference of order Nfr /LP' between the value of the integrated response function 
in the thermodynamic limit and the value at finite size. 

We can also look at the link with the static analysis. Like in the case of 
mean field spin glasses we can suppose, following [D^, that for i, t^, — > oo, 
C{t,tw) q and X[C{t,tio)] x{q), where x{q) is the cumulative distribution 
oiPiq): 

x{q) = r dq' P{q') (12) 
JO 

If this link is valid we can connect the susceptibility multiplied by the heat-bath 
temperature at a certain correlation value {S[C]) to the integral J^dqx(q). 

In our case there is an agreement between dynamic data and the values of this 
last integral computed on the static data. The two approaches are consistent 
if the external field is not too large, in order to avoid the non-linear effects 
that are neglected in the derivation of formula ^ The field cannot even be too 
small, though, in order to make the system move from the meta-stable state in 
which it has gone during tw The lower is the temperature at which we cool 
the system, the bigger is the external field necessary to make it explore other 
parts of the space of states. A probe at too low temperatures then is not really 
working because non-linear effects interfere heavily or because the system does 
not reach the aging regime. 
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Figure 10: Comparison between S[C] computed on dynamic and static data. The 
static curve is computed from the data of a system of linear size L = 32 at equilibrium. 
The other curves from the data of a system with 128 x 128 tiles out of equilibrium. 
Measurements starting at different waiting times tw are shown. The two plots above 
represent the situation at T = 0.35, where a perturbative external field of magnitude 
ho = 0.2 has been applied at t^. The two below show the relation between response 
and correlation functions at T = 0.30, with ho = 0.1. For large tw the dynamically 
determined S[C] tend to lay on the asymptotic static curve. On the rightside plots of 
the detail of the beginning of the bending is shown. 



Two examples of the agreement between dynamic and static data are shown 
in figure |l^. 
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4 Conclusions 



We have studied the thermodynamics of a tihng model built by Wang tiles. 
For this kind of system we have found evidence of a phase transition from a 
completely disordered phase, in which the tiles on the plane are completely 
uncorrelated between each other, to a phase in which they begin to present an 
organized, also if very complicated, structure. For T this structure becomes 
an exactly matched tiling. 

In order to characterize the phase of the system we have determined an 
order parameter with a non-zero mean value below Tc- It is the overlap between 
two tilings built with the same boundary conditions on the lower and the left 
sides of the lattices. Our data hint that it could have a non-trivial probability 
distribution under the phase transition. 

Under the critical point the tiling system shows the aging phenomenon: the 
answer of the system to an external perturbation and the time autocorrelation 
function depend on the history of the system. This brings to a violation of the 
fluctuation-dissipation theorem for T < Tc. The integrated response function 
times the temperature {S[C]) shows a progressive bending when the system 
leaves the stationary regime until it reaches a constant value {X = 0) and 
the dynamically and statically determined q values increase continuously with 
decreasing temperature. 

From this kind of behaviour it is not clear wether the model belongs to 
the class of systems showing domain growth or it is rather more similar to a 
spin glass in magnetic field. Indeed for very long times (small values of the 
correlation function) the fluctuation-dissipation ratio goes eventually to zero 
and it can not be excluded that the dynamics evolves through domains growth 
[ p^ , even though in our case the nature of the domains of tiles should still be 
theoretically understood. Nevertheless for a very large interval of time, i.e. of 
values of C in the plot of figur e {wj , the S{C) is continuosly bending (0 < X < 1) 
like in a spin glass model [p| |l l| , especially in the regions 0.3 < C < 0.6, for 
T = 0.35, and 0.4 < C < 0.7, for T = 0.3, as shown in flgure For values 
of C smaller than these the S{C) curves flatten but the predictions obtained 
from the static behaviour (the P{q) are shown in figure ||) are also nearly flat, 
making quite difficult the distinction between a domain growth dynamics, where 
the response function is constant in the aging regime, and a more complicated 
behaviour where even the response function shows aging. 

Eventually we show that our data are consistent with the equivalence be- 
tween the equilibrium function dq x{q) and the dynamically determined S[C] 
as t, tw oo. 
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